Assignment03

学生：麦泽霖 (12132211)

批改人：Juan-Li

得分： 49/50

PS3

1. Niño 3.4 index

import numpy as np
import pandas as pd
import xarray as xr
from matplotlib import pyplot as plt
%matplotlib inline

D:\Programs\anaconda3\lib\site-packages\xarray\backends\cfgrib_.py:27: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message
  warnings.warn(

da = xr.open_dataset('NOAA_NCDC_ERSST_v3b_SST.nc',engine='netcdf4')
da

D:\Programs\anaconda3\lib\site-packages\xarray\backends\plugins.py:61: RuntimeWarning: Engine 'cfgrib' loading failed:
Cannot find the ecCodes library
  warnings.warn(f"Engine {name!r} loading failed:\n{ex}", RuntimeWarning)

<xarray.Dataset>
Dimensions:  (lat: 89, lon: 180, time: 684)
Coordinates:
  * lat      (lat) float32 -88.0 -86.0 -84.0 -82.0 -80.0 ... 82.0 84.0 86.0 88.0
  * lon      (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
  * time     (time) datetime64[ns] 1960-01-15 1960-02-15 ... 2016-12-15
Data variables:
    sst      (time, lat, lon) float32 ...
Attributes:
    Conventions:  IRIDL
    source:       https://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/...
    history:      extracted and cleaned by Ryan Abernathey for Research Compu...

1.1 [5 points] Compute monthly climatology for SST from Niño 3.4 region, and subtract climatology from SST time series to obtain anomalies.

weights = np.cos(np.deg2rad(da.lat))
sst_weighted = da.sst.sel(lat=slice(-5,5),lon=slice(190,240)).weighted(weights)
region_mean = sst_weighted.mean(dim=('lon','lat'))
monthly_sst = region_mean.groupby('time.month').mean()
x = region_mean.size
anoma = np.zeros(x-2)
for i in range(0,x-2):
    k = i % 12
    anoma[i] = region_mean[i] - monthly_sst[k]
    
##  The criteria, that is often used to classify El Niño episodes, 
# is that five consecutive 3-month running mean SST anomalies exceed the threshold.
def moving_average(x,w):
    return np.convolve(x, np.ones(w),'valid') / w
mean_anomalies = moving_average(anoma,3)
mean_anomalies

# create a dataframe
# 'JFM','FMA','MAM','AMJ','MJJ','JJA','JAS','OND','NDJ','DJF'
alist = ['JFM','FMA','MAM','AMJ','MJJ','JJA','JAS','OND','NDJ','DJF']
templist=[]
for i in range(0,68):
    for j in range(0,10):
        templist.append(alist[j])
templist.append('JFM')
templist.append('FMA')
df = pd.DataFrame(mean_anomalies,columns=['mean_3_anomalies'])
df['mo'] = pd.DataFrame(templist)
df['date'] = pd.DataFrame(da.time.values)
sub = np.zeros(680)
df

	mean_3_anomalies	mo	date
0	-0.352137	JFM	1960-01-15
1	-0.307922	FMA	1960-02-15
2	-0.210943	MAM	1960-03-15
3	-0.240803	AMJ	1960-04-15
4	-0.225803	MJJ	1960-05-15
...	...	...	...
675	0.500900	JJA	2016-04-15
676	-0.072014	JAS	2016-05-15
677	-0.442696	OND	2016-06-15
678	-0.618626	NDJ	2016-07-15
679	-0.728378	DJF	2016-08-15

680 rows × 3 columns

1.2 [5 points] Visualize the computed Niño 3.4.

# Visualize the computed Niño 3.4.
df.plot(x="date",y="mean_3_anomalies")
plt.axhline(y=0.5,ls="-",c="red")#添加水平直线
plt.axhline(y=0,ls=":",c="black")#添加水平直线
plt.axhline(y=-0.5,ls="-",c="blue")#添加水平直线
plt.xticks(rotation=45)
plt.ylabel('sst_anomalies[℃]')
plt.xlabel('year')
plt.title('Niño 3.4(5°N-5°S,120°W-170°W)')
plt.show()

2. Earth’s energy budget

d = xr.open_dataset('CERES_EBAF-TOA_200003-201701.nc',engine='netcdf4')

2.1 [5 points] Make a 2D plot of the time-mean TOA longwave, shortwave, and solar radiation for all-sky conditions. Add up the three variables above and verify (visually) that they are equivalent to the TOA net flux.

fig = plt.figure(figsize=(10,6), dpi=120)
grid = plt.GridSpec(5, 9)  # 3 rows 3 cols

plt.subplot(grid[0:2, 0:3])
ax1 = d.toa_sw_all_mon.mean(dim='time').plot()
plt.title('TOA Shortwave Flux - All-Sky',fontsize=10)

plt.subplot(grid[0:2, 3:6])
ax2 = d.toa_lw_all_mon.mean(dim='time').plot()
plt.title('TOA Longwave Flux - All-Sky',fontsize=10)

plt.subplot(grid[0:2, 6:9])
ax3 = d.solar_mon.mean(dim='time').plot()
plt.title('Incoming Solar Flux',fontsize=10)

plt.subplot(grid[2:5, 0:5])
ax4 = d.toa_net_all_mon.mean(dim='time').plot()
plt.title('TOA Net Flux - All-Sky')

plt.subplot(grid[2:5, 4:10])
ax5 = (d.solar_mon.mean(dim='time') - d.toa_sw_all_mon.mean(dim='time') - d.toa_lw_all_mon.mean(dim='time') ).plot()
plt.title(' Add up three variables')

plt.tight_layout()

2.2 [10 points] Calculate and verify that the TOA incoming solar, outgoing longwave, and outgoing shortwave approximately match up with the cartoon above.

[Hint: Consider calculating the area of each grid]

temp=[]
x = d.solar_mon.mean(dim='time') - d.toa_sw_all_mon.mean(dim='time') - d.toa_lw_all_mon.mean(dim='time')
y = d.toa_net_all_mon.mean(dim='time')
a = plt.scatter(x,y)
plt.xlabel('calculate_toa',fontsize=15)
plt.ylabel('toa_net',fontsize=15)
(x-y).sum()

<xarray.DataArray ()>
array(1.6990367, dtype=float32)

incoming_solar:340

outgoing_longwave: 240

outgoing_shortwave: 99

2.3 [5 points] Calculate and plot the total amount of net radiation in each 1-degree latitude band. Label with correct units.

weights = np.cos(np.deg2rad(d.lat))
net_weighted = d.toa_net_all_mon.weighted(weights)
region_mean = net_weighted.sum(dim=['lon','time']).plot()
plt.xlabel('lat/°',fontsize=15)
plt.ylabel('toa_net_all_mon/ W',fontsize=15)
plt.title('Total amount of net radiation')
plt.show()

2.4 [5 points] Calculate and plot composites of time-mean outgoing shortwave and longwave radiation for low and high cloud area regions. Here we define low cloud area as ≤25% and high cloud area as ≥75%. Your results should be 2D maps.

fig = plt.figure(figsize=(10,6), dpi=120)
grid = plt.GridSpec(4, 4)  # 4 rows 4 cols
plt.subplot(grid[0:2, 0:2])
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon<=25)).mean(dim='time').plot()
plt.title('sw-low cloud area regions')
plt.subplot(grid[0:2, 2:4])
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon>=75)).mean(dim='time').plot()
plt.title('sw-high cloud area regions')
plt.subplot(grid[2:4, 0:2])
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon<=25)).mean(dim='time').plot()
plt.title('lw-low cloud area regions')
plt.subplot(grid[2:4, 2:4])
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon>=75)).mean(dim='time').plot()
plt.title('lw-high cloud area regions')

plt.tight_layout()

2.5 [5 points] Calculate the global mean values of shortwave and longwave radiation, composited in high and low cloud regions. What is the overall effect of clouds on shortwave and longwave radiation?

fig = plt.figure(figsize=(18,6), dpi=600)
grid = plt.GridSpec(1, 2)  # 1 rows 2 cols
plt.subplot(grid[0, 0])
d.toa_sw_all_mon.weighted(weights).mean(dim=['lon','lat']).plot(label='global')
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon>=75)).weighted(weights).mean(dim=['lon','lat']).plot(label='high cloud region')
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon<=25)).weighted(weights).mean(dim=['lon','lat']).plot(label='low cloud region')
plt.legend(loc='best')
plt.xlabel('time',fontsize=18)
plt.ylabel('toa_sw_all_mon',fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.subplot(grid[0, 1])
d.toa_lw_all_mon.weighted(weights).mean(dim=['lon','lat']).plot(label='global')
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon>=75)).weighted(weights).mean(dim=['lon','lat']).plot(label='high cloud region')
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon<=25)).weighted(weights).mean(dim=['lon','lat']).plot(label='low cloud region')
plt.legend(loc='best')
plt.xlabel('time',fontsize=18)
plt.ylabel('toa_lw_all_mon',fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.tight_layout()

3. Explore a netCDF dataset

T = xr.open_dataset('2016-2020-T.nc',engine="netcdf4")
T.t2m.mean(dim=('latitude', 'longitude'))

<xarray.DataArray 't2m' (time: 43848)>
array([286.67023, 288.65994, 289.85648, ..., 278.07663, 277.76575,
       277.50546], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 2016-01-01 ... 2020-12-31T23:00:00

3.1 [5 points] Plot a time series of a certain variable with monthly seasonal cycle removed.

weight = np.cos(np.deg2rad(T.latitude))
T_weighted = T.t2m.weighted(weight)
region_mean = T_weighted.mean(dim=('longitude','latitude'))
a = region_mean.values
index = pd.date_range('2016-01-01', periods=1827, freq='D')
k = 0
t_day = np.zeros(1827)
for i in range(0,1827):
    t_day[i]=a[k:k+24].mean()
    k+=24
tempdf = pd.DataFrame(t_day,columns=['t2m'])
tempdf['date'] = index
tempdf['year'] = pd.DatetimeIndex(tempdf.date).year
tempdf['month'] = pd.DatetimeIndex(tempdf.date).month
t_month = tempdf.groupby(['year','month']).mean()


monthly_T = region_mean.groupby('time.month').mean()
sea = np.zeros(60)
for i in range(0,5):
    sea[12*i:12*i+12] = monthly_T

certain = t_month['t2m'] - sea
certain.plot()
plt.ylabel('temperature[℃]')

Text(0, 0.5, 'temperature[℃]')

3.2 [5 points] Make at least 5 different plots using the dataset.

T.t2m.mean(dim='time').plot()

<matplotlib.collections.QuadMesh at 0x1fcf932bc10>

import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(0)

# example data
mu = 100  # mean of distribution
sigma = 15  # standard deviation of distribution
x = a

num_bins = 50

fig, ax = plt.subplots()

# the histogram of the data
n, bins, patches = ax.hist(x, num_bins, density=1)

# Tweak spacing to prevent clipping of ylabel
plt.title('Histogram of T: $\mu=100$, $\sigma=15$')
plt.xlabel('temperature')
plt.ylabel('probability dendity')
fig.tight_layout()
plt.show()

data = a
# basic plot
plt.boxplot(data)

{'whiskers': [<matplotlib.lines.Line2D at 0x1fc8d7f9610>,
  <matplotlib.lines.Line2D at 0x1fc8d7f99d0>],
 'caps': [<matplotlib.lines.Line2D at 0x1fc8d7f9d60>,
  <matplotlib.lines.Line2D at 0x1fc8d807130>],
 'boxes': [<matplotlib.lines.Line2D at 0x1fc8d7f93d0>],
 'medians': [<matplotlib.lines.Line2D at 0x1fc8d8074c0>],
 'fliers': [<matplotlib.lines.Line2D at 0x1fc8d807850>],
 'means': []}

T_weighted.mean(dim=('longitude','latitude')).plot()
plt.ylabel('t2m[K]')

Text(0, 0.5, 't2m[K]')

monthly_T.plot()
plt.ylabel('t2m[K]')

Text(0, 0.5, 't2m[K]')